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Abstract 


iUel depletion problem in a nuclear reactor is a non ! in ear 
constrained initial value problem. Ihe present study makes an 
analytic approach for the solution of depletion problem for a 
slab reactor. The methods adopted here are called seni— analy- 
tical in the sense that although closed-fora solutions are 
obtained in the form of series, different coefficients occur- 
ing in the series are to be evaluated numerically. Snall-tiue- 
step solution provides a rough but quick estimate of nuclear 
parameters - which is veiy important for preliminary nuclear 
power plant design. In the variational formulation with spatial 
power smoothening as desirable criteria, it is shorn how the 
performance index may be constructed and simplified to avoid 
nonlinearity. A technique for handling non-linear performance 
index is found out in the minimal control effort formulation of 
the depletion problem. Galerkin Method of weighted residuals is 
used to find the relation between the coefficients of expansions 
of trial functions for the state variables. Study of the deple- 
tion problem after refuelling is also done. A method has been 
indicated to find out the refuelling tine for achieving maximum 
benefit. Finally, some recommendations are made for the impro- 
vement of the problem and its solution. 
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1.1 in important item in the cost of fuel per kilowatt hour of 
electricity produced by a nuclear power plant is the fuel cost 
namely expenses incurred in the processing, the fabrication of 
new fuel elemaats and refuelling. Except in the case of reactors 
like CANTO where continuous fuelling is possible, reactors are 
normally shut down for refuelling resulting in less utilization. 

So from economic point of view, bumup of fuel should be as large 
as permissible with proper limitations set by nuclear parameters 
namely reactivity and also from mechanical and metallurgical 
considerations. The effect of irradiation puts a limiting factor 
to the bumup by causing dimensional change and physical damage to 

both fuel materials and cladding. When bumup is not restricted by 
radiation damage, the net reactivity change due to the consumption 
of fissile material (bumup), the accumulation of fission products 
and the buildup of isotopes, the capture of neutrons by control 
materials and other poisons determines the life time of the fuel 
elements. From the above discussion bumup (fissile material 
depletion) study is very important for an economic nuclear power 
plant desigi. 

(l) 

1.2 Martin Becker solved the depletion problem in a slab 
nuclear reactor by Least Square Variational technique and showed 
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■that his method, gave better results compared to other variational 

methods namely Galerkin method, collocation method, Ritz method. 

(2) 

Dougherty and Shen applied a semi-direct variational method to 
determine the time dependent coefficients of a modal expansion 

fb ) 

of the neutron fluxes for the multigroup kinetic equations. Kyong 
has presented a direct technique for the solution of a class of 
optimal control problems involving a distributed parameter reactor 
and a new method for solving simultaneous Fredholm integral equation 
of the second kind which arisesin the optimal control problem. 

(4) 

ICoyose and Suzuki have shov.a by using optimization methods in the 

fuel management and bumup studies of nuclear power reactors, that 

the requirement for maximum fuel bumup does not lead to power 

(5 ) 

distribution flattening. In the study of Fuller a trial solution 

is formed for the neutron flux by making expansions in known spatially 

dependent functions called trial functions. The undetermined time- 

dependent functions, called amplitude functions are then found by 

using the weighted residual procedure. Kessler^ ^ has used the 

time discontinuous synthesis method to describe the space-dependent 

dynamics behaviour of fast reactors. The space dependency of the 

neutron flux, the material temperature and the feedbacks are 

treated in cylindrical geometry. Two dimensional trial functions 

for the multigroup fluxes and the concentrations of the precursors 

are determined in an iterative procedure. The methods of variational 

( 7 ) 

synthesis are applied by Stacey to the problem of computing the 
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optimal control for spatially dependent reactor models, furthermore 
a generalized formalism is developed using one group neutron diffusion 

( s') 

theory, with and without delayed neutrons. Temey and j?cnech v 

applied dynamic progra mmin g and a direct flux synthesis method to 

determine the optimal control rod programming which minimizes power 

peaking throughout core life for a pressurized-water reactor (F II?.) 

.(9) 

Motoda and Kawai discussed a three dimaisional multigroup diffusion 

theory code for use in fast reactor analysis. The code can be useo 

for computation of material bumup and fission product buildup for 

specified time intervals. A variational treatment of the bumup 

optimization of continuous refuelling is presented and numerical 

solutions are given for a slab reactor by Motoda^ 10 ^. Kalimullah^ "* 1 

applied synthesis programming to solve some problems of optimal 

control-rod programming of a nuclear reactor discribed by the point 

kinetics model, furthermore, he presented four computer codes to 

( 12 ) 

solve different kinetics problems by synthesis method, Mate Ison 

discussed in his paper two schemes in which solutions to the three 

dimensional transport equation can be synthesized from two dimensio. 1 

transport solutions. Derivations are presented which employ a 

weighted residual technique applied to the second order form of 

(l3) 

the transport equation. 3h Stacey 7 s paper the developments of 
variational principles that admit' discontinuous trial functions 
which need satisfy neither the final and initial conditions nor 
the external boundary conditions of the physical prohLem are reviewed 



the mathematxe- 


and generalized. He discussed in his later paper 
difficulties that arise when discontinuous trial functions are 
substituted for continuous functions. These discontinuous trial 
functions are investigated by formulating the p rob Ian in terms of 
wellknown step functions and their derivatives. 

1.3 3y studying different works till now the least Square 
Variational method is found to be a powerful technique for solving 
reactor problems. But this method in general leads to nonlinear 
equations which are time consuming to solve. Hence in the present 
work suitable construction of the performance index is made leading 
to a set of linear algebraic equations. Most of the work in the 
fuel deplection problem uses maximum bumup as the main objective 
while attempts are made here with new objectives as spatial power 
smoothening and minimal control effort. The nonlinearity in the 
performance ind§x is removed either by proper substitution or the 
use of linearisation technique of the actual nonlinear problem. The 
construction of trial function is made in such a way that several 
functions involved in the process are analytically integrable using 
orthogonal property. Thus numerical integration is avoided and 
hence computation time is saved. 

1.4 The fuel-depletion problem in a slab reactor of thickness 2a 
is considered for the present study. The time-dependent one— group 
diffusion equation is given by, 



{D(x,t) $(x,t)} + (vo f - o a )u( x,t) <}(x,t) - I ap (x,-fc) cf (x,t) 

“vft (x ’ t > t 1 - 1 ' 

where , 

= thermal neutron flux 
= diffusion coefficient 
= average number of neutrons per fission 
= fission crosssection 
= absorption crosssection 
= atom density 

= macroscopic absorption crosssection of poison 
material 

= speed of thermal neutron. 

The arguments will be omitted in next equations for 

In fuel-depletion problem, the significant change in 
flux and power distribution takes place over the course of several 
months. So, the time derivative in equation (l.l) can be safely 
omitted since flux variations in small time intervals are negligible- 
Thus, for practical consideration the above equation (l.l) may b- 
treated as a steady state equation, 


b(x,t) 

i(x,t) 

V 

a f 

a 

a 

N(x,t) 

1 (x,t) 

ap 

V 

simplicity. 



3b. 


) + ( vcr f 


0 

a 


) N<t - I 


ap 


b = 0 


(l-2) 



The variation of fissile isotope concentration is repre- 


sented. by the following equation 


3N 

at 


a t N 
a 


(1.3) 


In the present problem it is assumed that there is only one 
fissile isotope and there is no buildup of fission products. -tLsOr 
one group diffusion theoiy is assumed. Initial flux and fuel concen- 
tration are given at the start of the reactor. The depletion problem 
will be solved with the diffusion theoiy boundary condition that the 
neutron flux goes to zero on the boundary of the core. 


finally the constraint of the problem is that the reactor 
operates at a specified power I ’ (t ) . Mathematically, it is given as 
an integral constraint involving only the space variable, 
a 

G o f / dx N $ = P'(-fc) 
o 

i.e., / dx K $ = P(t) 0*4) 

o 

where G = energy produced per fission 
P(t) = P f (t)/& ° f 

P(t) is given as a polynomial of degree M— 1 

P(t) = ? P^' 1 0-5) 

m=1 

Generally reactors are designed for constant power operation > 
but variations upto a maximum of 10/2 can be allowed. Beyond the 1C£: 
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■variation a reactor is normally shut down for safely. Thus if is 

the initial reactor power then the contribution of other terms of 

the polynomial never exceed 10$ of P . 

1 

Summarizing, it can be said that the fuel depletion problem 
is a nonlinear constrained initial value problem in the form of 
partial differential equations with an integral constraint. 

1.5 Chapter 2 makes an attempt for the solution of the depletion 
problem in small time steps by the assumption of neutron flux remain- 
ing constant over each small time step. Due to this assumption the 
solution is greatly simplified both derivationwise and computation- 
wise. It will be shown later that the results obtained are fairly 
accurate compared to existing methods. Chapters 3,4,5 are variational 

formulations of the depletion problem. The present problem is similar 

(15) 

to that taken by Martin Becker in his doctoral thesis . In Becker : ? 
least square variational formulation solution of a system of non- 
linear algebraic equations is involved, whereas in the optimal 
control method adopted here the functional can be so constructed 
that the problem reduces to the solution of a set of linear algebraic 
equations. Thus the method adopted here is better from the compu- 
tational point of view. The finite element method may also be 
adopted to solve the depletion problem. Although the desired accuracy 
may be achieved, computer time and data management may cause serious 
problem in the EEM. Another difficulty is that different approximations 
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in the PEM should be tested to know exactly which of these approidnrti „s 
gives the most accurate result with reasonable computer time. 

Chapters 3 >4 deal with an optimal control problem desired 
criteria being spatial power smoothening. 3h Chapter 4 fuel refuelling 
studies are made. Phe effect of refuelling decisions on the operator g 
tine and hence on the fuel cost is also given in Chapter 4. Chapter 5 

deals with a minimal control effort prOblem-f or ^constant power output. 

A technique for handling nonlinear perfonaance index is indicated in 
this chapter. 



CHASTER - 2 


In this chapter, the depletion problem is solved numerically 
in small time steps assuming that neutron flux in each time step 
remains constant. The governing differential equations and 
constraint (constant power) equation are rewritten here for easy 
reference. 


-2, 

D — g + a f “ ° a ) 

ax 


9N xkt 

3t =' "a*” 


I ♦ = 0 

L ap 


( 2 . 1 ) 


(2.2; 


/ dx = P(t) 

o 


(2.3) 


Symbols have their usual meaning as described in chapter 1. These 
equations are solved with the following boundary conditions 


and 


<j>(x,0) = 6 (x) 


N(x,0) = N 


A n TTX 

A Cos — — t 

o 2a 


(2.4) 


where A and N are constants, 
o 0 


Also, neutron flux is zero on the boundary of the core i.e., 


<j>(a,t) = 0 V t 


(2.5) 


This kind of small time-step-solution is a very common and 
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simplified approach to solve nonlinear partial differential equations. 

Although it is not a new technique, such an explicit solution for 

flux and atom density has not been sofar applied in the study of 

( 16 ) 

fuel depletion problems* The method was suggested by Zweifel 
in his text book on Reactor Physics for the bumup problem. 

If the reactor operates at constant flux, <j> is independent 
of time* However, reactors generally operate at constant power, 

since N changes with time, then $ will also change with time* 

That is, <j> is really a function of N, because as the uranium bums 
the resultant decrease in macroscopic cross-section will tend to 
make the flux rise to keep the product constant* 

The long-term time-dependent reactor equations are Integrated 
in small time steps, in each of which $ is assumed independent of t. 
Then at each step a new N is calculated, the corresponding $ is 
calculated and the procedure is repeated* One difficulty is that 
as the fuel is depleted, the control rods will be withdrawn to 
mai n tain a critical condition* This tends to change the flux 
distribution in the reactor. However, the short-time-step approxi 
matron is generally adequate to account for such slow variation of 
flux with time* 

2*2 The fissile number density as given by equation 2.2 is, 

~ = - <x ■ * at 
IT a 



integration and using the initial condition, 

t 1 

- / <? a <J> (x,t) dt (2.6) 

uCxjt^ = II e 

If <fr(x,t ) is assumed to be constant and equal to its initial 
flux value for the entire first time step t^= At after the 
start of the reactor, then U(x, At) is given by, 

, , - c a * 0 « At 

U(x, At) =N o e (2.7) 

Then from the constant power constraint equation 2.3, 
a 

/ dx U(x, At) <j»(x, At) = E( At) (2.8) 

o 

Substituting the spatial variation of <f>(x, At) by a cosine 
function, which is the solution of the neutron diffusion equation 
in a bare reactor and the time dependence by an amplitude function 
the equation 2.8 becomes 

a 

/ dx A ( 0 +At) C0S 2a At ^ = *(**) (2*9) 

where ♦(*, it) = A (q ^ 4t) ooe g 

The amplitude ^ ^ ^ will be a constant over the time 

stop At, but changes in the subsequent time intervals. 
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For the first time step from the above equation, the amplitude 
function is given by, 

l(At) 


l (0 -v At) ““ a 

/ dx N(x, At) cos ^ 
o 2a 


( 2 . 10 ) 


The integration in the denominator is carried-out numerically by 
dividing the reactor into well defined zones and them applying 
Simpson’s Buie. Thus the amplitude term ^ ^ ^ is evaluated 

for the first time interval* 


The flux 4>(x, At) thus evaluated, is assumed constant for 
the next interval between At and 2At. Fissile atom density 
Iff (x, 2At) is given by, 

a b( x > At) . At 


N(x, 2 At) = N(x, At) e 


a T 


( 2 . 11 , 


By repeating the procedure described above, the amplitude term 
A(^t ->-2At) 116)106 2At) can be evaluated. 

The poison cross-section at each time step is calculated 
from the neutron balance equation 2.1 as, 


■ J ap = fc Cl> ]/ * + (v ° t ~ H 


J f a' 


( 2 . 12 ) 


The procedure is repeated to calculate <}> and AT at each time 

step until £ is reduced to a certain preset value, 
ap 


When 
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the preset value is reached the reactor cannot operate any longer 

and hence refuelling becomes necessaiy. Since the present work 

( 15 ) 

was compared with Becker's Least Square variational method 
where the preset value of £ was not available* the problem was 
continued for various time intervals upto a maximum of 300 days. 
The final bumup data was obtained at all time steps upto the 
maximum of 300 days. 


2.3 Results of the depletion problem for a slab reactor of width 

200 cm fueled with U __ c ( = .4 ) with initial multiplication 

235 *28 

factor of 1.032 are discussed in this section. To compare the 

solutions of the present work with existing solutions available 

(l5 ) 

in literature, the Input parameters used by Becker have been 
taken for the present problem. 


D = 1 cm (constant) 


*■*24 2 

a = 200 bams (l bam =10 cm ) 
a 


va f = 400 bams 

3Q-4—I 

l(t) = 6.9 x 10 cm sec 


= A 0 cos g 


13 -2 -1 

where, A = 6.284 x 10 cm sec 
7 o 

19 -3 

IT = 2 x 10 cm 
o 

l +_n = 3.7605 10" 3 cm -1 

%p t=0 


At = 0.5 day 



Etie table 2.1 shows the neutron flux distribution and 

fuel distribution oxer space for different reactor operating 

times for both the present work and Becker's work. These 

results are plotted in figure 2.1. The figure 2.2 shows the 

variation of poison cross-section with reactor operating times. 

Note that the poison cross-section goes down steadily because 

238 239 

of the assumption that there is no conversion of U to ?u 
235 

and only U is the fissile material which is constantly 
depleted. Becker's results are based on the principle of minimizing 
the least square error of accurate solution and the approximate 
trial function (synthesis) method. The diffusion theory boundary 
conditions are applied in the present work and as well as in Becker’s 
analysis. The finite difference results of this chapter are compared 
with Becker's work only to establish the credibility. Later on the 
poison cross section variation obtained in this chapter is compared 
with the cross section obtained based on minimal effort criterion. 
Although the method adopted here is simplified and approximate, 
the results are close to Becker's solution within a maximum error of 
13$. The discrepancy is due to the assumption of a cosine function 
for spatial flux variation as opposed to a trial function in the 
form of a series in Becker's problem and also due to the assumption 
that flux remains constant in each time step whereas flux is a 
continuous function of time in Becker's method. 
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CHAPTER - 3 


3.1 la this chapter an attempt is made to attack the depletion 
problem as an optimal control problem and the solution is achieved 
by the technique of calculus of variation. The equations concerned 
are stated again for easy reference. 

'4* Cvo f -o a m-I ro *.0 (3.1) 

oX 

( 3 - 2 ) 

a 

/ dx N<j> = P(t) (3.3) 

o 

These equations are to be solved with the same boundary conditions 
as in chapter 2. 

In the language of optimal control theory, N and $ are called 

state variables and T is called the control variable. The 

^ap 

objective of the optimal control problem is to achieve smooth power 
density throughout the space over the entire operating time of the 
reactor. It is desirable to have smooth power density from the 
engineering design aspect since the core temperature will be uniform 
and constant and the total power output will be higher. The quadratic 
performance index for the present problem can be written mathematically 



where a a = weight factors, 

IT',’’ — power density. 

P(t)/a = average power density at time t. 

The weighting factors a, end ct^ are arbitrarily chosen to give the 
appropriate importance either to the variations in the power density 
or the gradient of power density. This quadratic performance index 
is to be minimized by proper selection of the state -variables N and y 

Note that the control variable I does not occur in 3.4 

ap 

explicitly . The control, variable appears only in the governing 

differential equation 3.1 which is came as diffusion equation. Since 

the control variable does not appear in the performance index, 

according to optimal control ':heoiy the performance index can be 

minimized satisfying the dynamics equation 3*2, the integral constrain 

3.4 and the initial and boundary conditions on the state variables. 

Thus by proper choice of the performance index the control variable 

may bo omitted from the main optimal control problem which will mate 

computation easier. When, the reduced problem is solved for N and .1 

by the calculus of variation oechnioue then the control variable 7 

L ao 

in found out easily from the neutron balance equation 3.1. 

3.2 The trial function which satisfies the initial and boundary 
condition for the state variable 4 is chosen as, 



20 


M 


<!> (x,t) = 4> ( x »^)l t _ 0 + I l a« cos-^ — t m 


t=0 - “ x, m 

a =1 nt=1 


2a 


L M 

= * 0 M + l i 

fi . =1 m =1 


kn 7T x 




m 


cos 


,1 


2a 


(3.5) 


a £m are expansion coefficients which are constant. 

In order to satisfy the boundaiy condition that the flux goes to zero 
at the outside surface for all t, k^ have to be odd numbers. Thus, 


k £ = 2Z - 1 ; SL= 1,2,3,... . 

This trial function also satisfies the initial condition that, 


<f> (x,0) = d o (x) = A q cos g . 

In the present problem M is chosen equal to 2 i.e. a second degree 
polynomial in time is taken for ease of computation. Choice of 1 
was based on less than 1% variation which led to 12 term expansion. 
Similarly the trial function for the state variable N is constructed 
as, 

HU,t)=N(x,t)| t=0+ l l 

£=1 m=1 


= et 


i 

■ l 

a =1 


M 

I 

m=1 


Jim 


cos 


k a ir: 

2a 


ja 


(3.6) 


where, b^ m = constant expansion coefficients 


and k^ = 21 - 1 j A = 1,2,3, 
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This trial function satisfies the Initial ccndition, 

N(x, 0 ) = IT . 

3.3 a) The performance index as given in equation 3.4 is very 
difficult to handle because of the nonlinear product term of the 
state variables N and <$ and hence product of the expansion coeff- 
icients a £m and b^ m will appear in the performance index. 


Tfom equation 3*2 it can be written 



1 ? N 




(3-T) 


1 

where, k = ~ . 

a 

Putting 3.7 in equation 3.4, the performance index will be dependent 
only on one state variable namely the number density IT. 


P.I. 


X CL 

/ / ‘ 




at 


P\2 




a ^ k h **■ dtJ 


(3.8) 


This performance index is easier to tackle since only one of the 
unknovyn coefficients of expansion namely appear in it* 

3 * 3 . b) The constant power density constraint is 
a 

/ dx N as p(t). 
o 

Once again eliminating the nonlinear product term J3L$ by using 3*7 
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- / dx (“|t) = *“(t) • 
o °a 

Using the polynomial expansion given in 1.5 for P(t) and the trial 
function for N as given in 3.6, 


JL 

a 


a 


/ 


a o 


Is 

l 

£ =1 


M 

I 

m=1 




cos 


k^TTX 

2a 


m t' 


m-1 


dx 


M 

- I 

m=l 


P f 

m 


n-1 


Nov/ equating coefficients of like powers of t, 


a 




P Y b. / cos 

m o .i £m J 
a £=1 o 


dx 


2a 

m — 1 j2j3»**»>M 


After carrying out the integration in 3*9, 


P 


2am 


hi £ 


m 


===- y b , 

ire L £m k 

a £ =1 £ 


m = 1,2,3 >••• M 


with little algebraic manipulation, 
I / .%£ 


('ll 

b, = H + y b 1 ir L - m= 1,2,3,.*. M 
1m m %Z 2 k £ 


(3-9) 


(3.10 


2am 

Thus the integral constraint is converted into a set of equivalent 
algebraic equations. The equation 3.10 tells that out of L x Ivl 
variables, M variables (b ^ , m - 1,2,...,M) are constrained and the 
rest (l x M - M) are independent variables. 
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3*3 c) The trial function for IT given in 3.6 after some algebraic 
manipulation becomes, 


1 M M 

H = *0 + I l b <5 + I H 003 ~ t 1 

1=2 m=1 *- n 1“ mi, “ 


where. 


= cos . 
iim 2a 


K t, X / . si 

“ t “ + i=H- 003 ^ t 1 ” 
\ 2S 


So, 


3N 


1 M 


I 


at " X I. V G , mt + I H m COS II * mtmH 


o amt 1 a 

X/=^ m=l m=l 


(3.r) 


and., 


3 /31T- 


I U 


M 


. v (If) = I y b a GL + + y H (- — )sin XZ mt m_1 
3x 3t £ £ 2 m =! ^ 2a ; 2a 


(3.12) 


G-^mt represents differentiation of G with respect to t 
and represents differentiation of G^ m with respect to t and 
x both. They are given as 


vx .j (-i ) itx m-i 

W 003 ir~ *Tr~ ccs -£ mt 


^Jlmtx ~ (*" 2a 


kg v x 


■) sin — mt 

' 2a 


-1 (-1 ) r it \ . irz j."**'- 1 

, J— ( — — ) 33JQ ^mt , 

k % v 2a' 2a 


After substitution of 3.11 and 3.12 in the P.I. as given by 
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I a 


5.1. 


ft p ^ M m 

-//{«! k ( ^ I b G + I H cos ^ mt a -1 + 1)2 
O o .£=2 m =1 2,111 tot ^ m 2 a + a' 


I M 


+ a o k 2 ( 1 l b 


M 


iU 2 m=1 2,111 &£mtx + E ^~ 2a ^ sjil 2a ) 2 > <^t 


... (3.13) 

lo obtain the optimal solution, partial derivatives of the 
performance index as given by 3.13 w.r.t. all the independent variables 
( b £ m » ^ = 2>3 >...jL $ n = 1 } 2 ,...,m) are equated to zero. 


B P* I. 

= 0 , £ ss 2,3, ... ,L $m = 1 ,2, • . . ,M • 

£m 


(3.14) 


Por the present problem the above conditions reduce to a set of 
( Ir-1 ) x M linear algebraic equations. 

After some algebraic simplification these equations take the form, 


k‘ 


“l b £m ^\at J G £mt^ + k “l ^ \ b.. (G , G ) 

1 M ±Z £mt* ijt y 

j ^ m 


o M 

+ "l k(P ’ & nmt } + “l k l < H m 003 g mt > W + 

mps f 


+ \ k 2 (G tmtT , + a, ll b <4 G, 44 J + 


'Unto:’ -imtat' ”2 " ,, "ij ' Mtx’ “ijti" 

X /= W 

j ^ m 


M 


+ <* V (h (— •— -*) sin mt 111 G ) — 0 

p L * m v Oa / ? U 0Tn+-%r / 


2 L N “m v 2a" 2a 
nt=1 


jimtx 

A* 3 2 ,3 9 * * * ,1> 
m = 1,2,. . .,1 


(3.15) 


k » .o-vL UBHaRY 


tm» Mo. 
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where, the inner product (x,y ) denotes integration of these variables 
over the entire space and time. It is interesting to note t ha t 
= 0 as proved below, 


la Ta koffx z 

/ / . mt 11 ' 1 + ML c00 g . hjO-I j 


o o 


0 0 


2a 


dx. dt 


JL). 

z 


The algebraic equations finally take the form 


°1 b £m^ G S,mt* + a 1 \ } b ij & ijt^ 


i ^ Z 
3 A 


M 


+ cx V (h cos ~~ mt E> ""* ,Gr nj _)+ a ( G- • G- ) 
1 m=1 m 2a 2 ^ Zntx 


+ « l l b (q , G. )+« y (h(-~) sin^mt 33 " 1 , 
2 xj £mtx T x^tx 2 m 2a 2a * 


j / a 


W'° ( 3 - 16 > 


& = 2,3 , • • • , Ii 
m = 1,2,. ..,M 


This set of linear algebraic equations are solved to obtain 
the independent variables b^. ( Z s 2,3,..., L ; m = 1,2,...,M). The 
unconstrained variables b^ m (m = 1,2,k.., M) may be calculated from the 
constraint relation 3.10. 



26 


When all b^ are known, the closed form solution for N(x,t) 
is given by equation 3.6. The neutron flux $(x,t) is calculated f nor- 
ths equation 


3N __ 

at = 


°a 45 


or, 


♦ - - <!!>/% * 


(3.17) 


!Ehe whole calculation is carried out for a certain reactor 
operating time T^ days. The poison cross-section, is found out when 
N and <t> are known from the equation 3.1 as 

2 

I + (v« -o )B . 

3x 

If the value of Z is greater than a present desigied value, the 
operating time is increased to days by an amount At. The procedure 
is repeated until the poison cross-section reduces to the preset desi^:. 
value. Since the end value of the poison cross-section was not 
available for the present problem, the reactor was operated upto a 
maxi mum of T = 325 days and the final bumup data was obtained at 
T = 325 days. 


3.4 Results of the depletion problem in a slab Nuclear Reactor 

N 

of width 200 cm fueled with U (rr*»- = .4) with initial multiplication 

^ ^28 

factor of 1. 032 obtained in this chapter are discussed in this section. 

( 15 ) 

Input parameters are same as in Becker’s problem 


only P is considered 
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as a first degree polynomial in time as opposed to a constant in 
Becker’s problem. 


D = 1 cm (constant) 



*«rv>er ° A =r 

* G 



200 bams (l bam = lo" 24 


400 bams 


A cos ~ 

O 2a 

6.284 >< 10 13 cm" 2 sec" 1 


2 x 10 ^ 


3.7605 x icf 3 cm"" 1 



T^ = 200 days 
At = 25 days 


P(t) = P. + P 2 t 

39 -a 

P.j = 6.9 x 10' can sec 
P 2 = P/(l0t) . 


The summation indices are N = 12 and M = 2. 

The weight factors are 0=1 and - 50 . 

The choice of the polynomial for power leads to a 10^ increase 
in power but this is done in order to facilitate the solution and 
to see how the nuclear parameters change with a polynomial power 
constraint. The table 3.1 shows the neutron flux distribution and 
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fuel distribution oyer the space for different operating tines and 
these are plotted in figure 3.1. The figure 3*2 shows the variation 
of poison cross-section with reactor operating tine T. The poison 
cross-section here also goes down steadily due to the assunption 
that no Ba 5 is produced during the operation of the reactor. The 
variation in N and is similar to the results of the previous 
chapter but magnitudes cannot be compared as the present case deals 
with spatial power smoothening throughout core life tine as opposed 
to constant flux in small tine step intervals of the previous 

chapter. There are no experimental results available for simple 
slab geometry investigated in this thesis * The input data was not 
varied to consider practical situations, due to lack cf computer 
availability . 
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FIG. 3-2 


CHAPTER - 4 


4 • 1 IUel management is essentially the study of fuel loading 
schedules. Some of these schemes may require changing the location 
of the fuel elements within the core according to a predetermined 
schedule. Although these schemes may increase the burnup and reduce 
the expenditure of fuel, they may add to the cost of reactor 
operation, especially if shutdown is necessary for rearranging the 
fuel. So, careful economic analysis has to be made to determine 
the merit of the various scheduling schemes proposed by optimization 
studies. In the present chapter general fuel loading schemes are 
described and batch loading scheme is investigated in studying the 
fuel depletion problem. 

Batch loading s 

The reactor is loaded uniformly with complete core at one 
time and this is irradiated without moving the fuel elements. When 
the system ceases to be critical the whole core is discharged and 
reloaded. As a modification of the simple batch irradiation scheme, 
the core may be divided into two or more radial zones, with fuel 
of different enrichment in each zone. 

Centre to Outside loading : 

Prom time to time, fresh fuel material is added to the 
centre of the core, where the leakage of neutrons is low and henco 
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their importance is high and progressively moved to outer radial 

positions. ®ie spent fuel is discharged from the outer radius of 
the core. 

Outside to Central loading : 

In this procedure, the fresh fuel is charged near the outer 
edge and progressively moved towards the centre from where it is dis- 
charged . 

Bidirectional loading : 

In bidirectional loading, charging and discharging of the 
fuel are achieved by pushing relatively- short fuel bundles through 
the core from one side of the reactor to the other, but from opposite 
directions in adjacent rows of fuel channels. 

4.2 She reactor is operated according to the desired criteria 
of spatial power smoothening as in chapter 3 in the present study. 

The same depletion problem is considered here also. The operatic a 
of the reactor is stopped after time T which varied from 200 days 
to 300 days in steps of 50 days for refuelling. The comparatively 
more depleted fuel zone ia replaced by fresh fuel of the uniform ato 1 - 
density corresponding to the concentration after T days of 
irradiation at the spatial location of x «* 60 cm as shown by the 


point B in figure 4.1 
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Olhe comparatively less depleted fuel zone (zone 2) is not 
changed. For further reactor operation the initial fuel configuration 
is now constant value in zone 1 as indicated by the IB line and the 
depleted values as given by the curve BC in zone 2. Mathematic ally 
this can be represented as, 

“ H(6O)l te2 00 * tD o« - °60 W> 

(after refuelling) + N(x)| t _ 20Q * U g0 (x) 4.1 

N sb New fissile-number density, 
where the Unit step function U is defined as 
U (x) = 0 x < 0 

U o (x) a 1 X > 0 
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v 60 W » o 

x < 60 

U 60 (x) ='1 

X > 60 


4.3 The new trial function for fuel density 


is given by 


NCx.t) . St*,0) ♦ l ib tnCOS !^ t m 
£=lw=l 2a 


4.2 


whore are the expansion coefficients. 

Poliowing the procedure outlined in chapter 3 the equivalent algebraic 
form of the integral constraint on power is given by 


b 1m = 


I f 

h + y b d4iU_ 
n Z=2 V 


SLm K. 


4.3 


m 




-P a ir 

where H = — — - 

m 2am 


Incorporating 4.3 in 4.2 and after some algebraic manipulation. 


^ _ L M M 

N(x,t) = N(x,0) + I l b G (x,t) + I H cos f 

£=2 m=1 q =1 n 2a 


where, 


^n' 000 2a 


K t* x n (-1)* ... 

t + — L — cos rr- t 


ttx ja 
2a 


4.4 


The spatial power smoothening being the desired criteria, 
the performance index remains the same as in chapter 3* Thus 


4.5 


P.I. = / / [ot (k 


3N 

9t 


+ £) 2 

a 


0>(K 


_3_ 

3x 


3^2 
9t 


) ] dx dt 


The function H(x,0) is continuous as seen in the figure 4.1, 
but its derivative ~ is not continuous. To avoid the point of 
discontinuity integrations are perfoimed upto a small neighborhood e 
of the point x= 60 cm. The performance index may then be written as 


p.i. = / T / 60 <«j(k h) 2 . at A ) 2 } 2 & dt 

0 0 A 3 t a i dx a t 


T 100 nN p 2 « gu 2 

* f L * a 2 tK lr c w' :)) tedt 

o 60+ e 1 3t a z 9 X 3t 4,6 


where e is a veiy small quantity and is taken to be 0.01 for th~ 
computations performed in this chapter. The performance index in this 
particular form does not produce any computational difficulty. 

The expansion coefficients b^ are solved by the same calculus 
of variation technique as in chapter 3- When b^ are solved, closed- * 
form solution for N is achieved through equation 4.2. The calculation 
of neutron flux and poison cross-section is exactly the sane as in 
chapter. 3. 

4.4 The reactor description and the input parameters are same as 
in chapter 3- 

The table 4.1 shows the neutron flux and fuel distribution 
over space for different operating times T after refuelling. 
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They are plotted in figure 4.2. The figure 4.3a shows the variation 
of poison cross-section with operating time T. is observed in earlier 
chapters, the poison cross-section steadily goes down with time. 

This is because of the fact that it is assumed that there is no 
plutonium production from U 

238 

The objective for refuelling may be two fold i) to Increase 
the bumup of the fuel ii) to increase the overall corelife of the reactor 
But to ascertain whether a particular fuel schedule should be adopted, 
careful economic analysis should be carried out first. In the present 
problem, the comparatively more depleted fuel-zone is replaced by a 
uniform, fuel density of higher concentration. The less depleted fuel 
zone near the core end is not replaced. So, this allows effective 
burning burning of the fuel near the end which is not realised if the 
whole core is refuelled. Due to the introduction of fresh fuel, the 
amount of poison cross-section for controlling criticality increases 
and hence the reactor may be operated in the critical condition for 
longer time. 


An approximate estimate of the fuel cost of the newly 
introduced fuel may be obtained in terms of reactor operation days. 
Prom fig. 4.1 it is given as 


Puel Cost 


Area enclosed by ABD 
Area enclosed by AEC 


x Time of refuelling 


4.7 


Por the present study, refuelling is done after 200, 250 , 300 days and 
the fuel costs of the newly introduced fuel in terms of days are 
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25, 31, 37.5 days respectively, How if the reactor is not refuelled 

the poison cross-section reduces to 1.901 * icf 3 on” 1 after 300 days 

of operation. As seen from figure 4,3, the reactor operates for 

(200 + 133) - 333 days in total when refuelling is done at T = 200 

days before it reaches the poison cross-section 1.901 x ■jq~^ 

Similarly, the reactor ■, operates for (250 + 88) = 358 days if 

refuelled at 1 = 250 days and (300 + 45 ) » 345 days if refuelled at 

T = 300 days before the poison cross— section goes down to 
-3 -1 

1.901 x 10 cu . Therefore, the reactor operates for 33, 38 a n d 45 
nays more if refuelled at T — 200 , 250, and 300 days respectively. 
The gain due to refuelling may be obtained by subtracting this 
extra operation days from the fuel cost of the new fuel introduced 
at the time of the refuelling. 

Thus the gain achieved due to the partial refuelling at 
T » 200 days is given as (33 - 25 ) = 8 days. Similarly the different 
gains for refuelling at T = 250, 300 days are 7 and 7.5 d^ys. 

The number of days gained vs. refuelling time may be plotted from 
where the optimum tine of refuelling can be chosen for maximum gain. 
In the present study the gain appears to be almost constant with 
time of refuelling according to Fig. 4.3- This is due to the rmgh 
estimate of fuel cost calculated on area basis and also due to the 
fact no plutonium production is considered; So in the present 
case it is very difficult to make a decision on refuelling tine for 
achieving a maximum gain of reactor operating time. 
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CHAPTER - 5 


5*1 This chapter makes an approach to the solution of the 
depletion problem as described mathematically by the following 
equations with a desired criteria of achieving minimal control 
effort* The minimal control effort is accomplished by variational 
technique as outlined in chapter 5* The governing equations are* 


g 2 a 

O’ ° f - V N* - l n ♦ . o 

(5.1) 

~ = - a AN 
at a Y 

(5.2) 

a 


/ dx N$ = P(t) 

(5.3) 


o 


The boundary conditions on TS and <f> are the same as in previous 
chapters * The performance index is chosen so as to minimize the 

parasitic capture of neutrons in the control material which will 
indirectly result in less fuel consumption. Mathematical this 
is given by, 


t a 


P -I- = / / lap 
o O 


dx dt 


(5.4) 
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where the poison cross-section i s the control variable. She 

fissile number density N and flux <j> are the state variables. 

The optimal control problem is to find the values of I and <l> 
for which the desired criteria given by equation 5.4 xs achieved. 
Mathematically, the p rob lan in this chapter is more difficult than tns 
problem posed in chapter 3 because the control variable appears 
explicitly in the performance index and hence nonline ar equation 
is to bo solved without further simplifications. Hie diffusion 
equation 5.1 may be used to find the desired values of fissile number 
density N and (J> by substituting these variables for in 

equation 5*4. Thus -Hie performance index changes to. 


P.I. 


T a D 2* 

j J + (v » o^)N} dx dt 

0 0 ^ 


(5,5) 


The nonlinearity arises because of the state variable <f> occuring 
in the denominator and this causes complexity in the present problei . 
Dote that the performance index in chapter 3 contained only the 
state variable N. But here it contains both the state variables 
N and 4 as given by 5.5* 

5.2 With the same basic philosophy as outlined in chapter 3, 
trial functions for the state variables f and H are constructed 
in this chapter.' They are-given'by equations 5.6 and 5.7. The 
constant power constraint as given by equation 5-3 puts restrictions 


*r4 


similar to equation 3-10 on the coefficients of expansion b^ . 
Following the steps as shown in chapter 3 this condition given by 
equation 5.8 is obtained. Here time variation of flux is chosen as 
linear as opposed to the quadratic polynomial chosen in chapter 3 
for mathematical simplicity . So, the governing equations are , 

Trial functions for and Ns 


<fr(x,t) = 


$ 0 M ♦ J 

1=1 


a cos 

a 


K 0 irx 

~ 1 

2a 


(5.6) 


I TTX 

N(x,t) = N o + l b^ cos t (5.7) 

1=1 

Equivalent algebraic equation for power constraint equation: 


b 


1 


P. TO 

_1 a 

2a 


+ 



iz ll 


1 K„ 


(5.8) 


Since both and N are involved in the performance index 5.5, it 

is essential to find a relation between the coefficients of expansion 

% 

a^ and bg_ in order to obtain complete solution. 


5.3 The Gale rk in method of weighted residuals is used to find out 
the relation between and b^ If the depletion equation 5.2 

is used without linearization in the Galerkin method, the resulting 
relation between a^ and b ^ will be nonlinear due to the product 
term, of the state variables N and (J> occuring in 5.2. To avoid 
this problem and to use the Galerkin method effectively linearization 
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of the depletion equation is done by assuming perturbation around 
initial values of the two state variables namely N o and <j> . Thus 

in the first order approximation the depletion equation may be 
written as 


3_(n 0 +*it) 

at + 


a a^o + ( N 0 + = 0 


where, dh a H - N 

c 

6<|> as tf> — <j) 


Neglecting the product of <5<J> and 6N , 


3N q 36N 

at*" + atT~ 


+ ° a 5N + c a N o + °a N o = 0 


(5.9) 

(5.10) 


Using 5.9 and 5.10 for 6 N and 5 <f> , the above equation reduces to 


3N 

at 


+ 


* N <f> +04 N 
a 0 a 0 


* N > 0 

a 0 o 


(5.11) 


The Galerkin method then yields the following equations, 



a 6 
a 0 



Kp ir x 

cos ) * 0 

2a 

P = 1,2,... ,L 


Kp DX 

where cos — are the basis vectors of the trial functions. 

The inner product (x,y) indicates integration over entire space 
and time. These inner products are easily calculated by employing 
orthogonality property of the cosine function. Except the first 
one they are not shown here. 
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(a a N o «, cos 


K wx T a 


— ) = / / {a N cos g. + ? a cos t] * 

oo ^ 2a 


V"\ 

-»v * 


X cos -i } dx dt 

^a 


2 CT a N o A o T 6 P i + a p i ~~ 


a 


where 


pi * Kronecker delta 


x.e. 


5 1 = 1 when p = 1 -j 


0 otherwise 


After carrying out the integrations and simplifying the algebraic 
steps, the final relation arrived at is, 


a Jlf !S " b £ (l +|) - 3a a I o A o 6^ (£=1,2,...l) (5.12) 

5.4 Using 5*12 and 5.8 in 5.6 and after some algebraic manipulation 
the expression for .flux becomes, 




<b + 


l “•>» S, + Y, 
£=2 


cos -g t - + A(f» 


(5.13) 


L 

where A<}> » ][ - bj, &£+ Y . cos ~ t (5.14j 

£=2 1 

The quantities Gj, and Y^ are defined as, 
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where suffix xx denotes double differentiation with respect to x. 

G fcra: is S iven as 
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Hence, 


d> d> + Ad) 

xx _ oxx xx 

<J> * 

4> 0 + Ad 

with the assumption Ad is sm al l compared to d Q (Ad << <J> 0 )» the 
expression in the right hand side can be approximated to, 
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and also using 5.14 and 5.15 with some algebraic rearrangements. 
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^xx ^ jj, Lb* 

~ \ seo <* 3 ^ l 2 T o G * ~ g> 


(5.16) 


where , 


r 2 = “ ^2af ~ Y l * @ 2 /A c 


Y = 1 - Yt/A 
3 1 ' o 


How from equations 5.7 and 5.8 after some algebraic simplification, 
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Prom 5.17 » 
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(5.18) 


where suffix t denotes differentiation w.r.t. t. 


Finally substituting 5.16 and 5.18 in the equation 5.5, the 
performance index now becomes, 

p.i. - / /V 2 + j - r G to seo {y 3 + j, r G * seo t' 
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(5.19) 


The solution of the minimal control problem is then given by, 


(a *2,3,. ..i) 


(5.20) 


This yields a set of linear algebraic equations since a linearized 
functional is used. These algebraic equations can be written as, 
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(5.21) 


i»5 j3 y « • « Ij 

The different integrals involved in 5.21 are carried -out numerically. 


The solution of these fc&-l) algebraic equations will give the 
coefficients of expansion bjj, (& =2,3,.**Ii). When b S( (S-=2,3,*..I<) are 
obtained b, la determined trom the constraint equation 5.8. Since 
all the expansion coefficients b £ are known, the closed form 
wolution for the state variable N is given by 5.7. solution 

for flux is obtained from 5.2 as, 


♦- - (!!>/>% 

The whole calculation is carried-oui for a reactor operating time 


T days 
5.1 as 


The poison cross-section is calculated from the equation 



D 


9 2 <f> 

ax 2 


+ (va - 


a a )n 


The operating time is increased every time by a certain n umb er of 
days and the whole calculation is repeated to calculate the poison 
cross-section for each operating time, until the poison cross-section 
reduces to a preset value. Since this preset value was not 
available, the reactor was operated for a maximum of 500 days and 
the final bumup data was obtained at T = 300 days. 


5.5 Results of the depletion problem as discussed in this chapter 

N 25 

for a slab reactor of width 200 cm flueld with — = *4) 

with initial multiplication factor of 1.032 are discussed in this 
section. 

Input parameters are same as in chapter 3* The only difference 
is that the trial functions in this chapter assume linear polynomial 
in time as against a second degree polynomial in chapter 3. The power 
is assumed constant as opposed to a polynomial in time in chapter 3* 
The table 5.1 shows the neutron flux and fuel density distribution 
over the space for different operating times T. They are plotted 
in fi gu re 5.1. The figure 5*2 shows the variation of poison cross- 
section with -reactor operating times T. As noticed in earlier 
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chapters, the poison cross-section steadily goes down with time. 

The reason is that TJ 5 is the only fissile material and it is 

assumed that there is no buildup of HJ 239 from U^q. Prom 

figure 5.2 it is seen that J ap is smaller in this case compared 

+o 7 obtained in chapter 2. It is therefore evident from 
L ap 

equation 5.1 that the fissile number density N is less compared to 
that obtained in chapter 2. Since N has gone down, the flux <J> 
has gone up to maintain a constant power as seen from equation 5.3* 

But since the present problem deals with an optimal control foxmulati.n, 
the results cannot be directly compared with the results of chapter 2. 


EL AND FLUX DISTRIBUTION 
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CHAPTER - 6 


T?ron the analysis carried-out in the present woik it is 
evident that the methods adopted here can tackle complicated non- 
linear initial or boundary value problem in the form of partial 
differential equations with constraint equation. 

But the depletion problem treated in the present study is 

greatly simplified. Some important factors such as conversion of 
238 239 

U J to Bu , build-up of fission products, presence of a reflector 
about the core, possibility of multiregion core are not considered 
here. So these refinements may be incorporated in the depletion 
problem for better understanding of the problem. 

In many nuclear reactors plutonium production at the begin- 
ning of life more then offsets u 235 consumption causing an initial 
increase in reactivity. The poison level should therefore increase 
to a maximum value before beginning to decrease. 

Taking plutonium production into consideration the governing 
equations for the depletion problem are changed as, 

Neutron balance equation : 

dLJ + (vcr f - ff a ) JH - I ap <J> + (v g o jg - a a9 ) V - 0 (6. 

3x 


where 




Mg = number density of plutonium 

o ss fission cross-section of plutonium 
fy 

Vg = average number of neutrons per fission of plutonium. 
Eynaraic equations : 


The rate equation for uranium is unchanged 




6.2 


The rate equation for plutonium is given by 


3M 

= - a <J) M_ + cr 0 N Q <!> 
3t a9 9 c8 8 


6.3 


where the first term is due to neutron absorption (fission plus 
capture) in Pu ovn and the second term is due to plutonium formation 

23 y 

from thermal neutron capture in U^g* 

°a9 = at,s0r P' fci0£1 cross-section pf Pu^g 
cr g a capture cross-section of TJ^g 
Ng = number density of U^g* 

As a first approximation, the concentration of TJ^g, which changes 
by a few per coat only during normal reactor operating time, nay 
be taken as constant# 


The constraint equation : 

a a .. 

/ dx M $ + a^g / dx Mg4> « p(t) 
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6.4 
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where, P(t) « p(t)/G. 

iill other symbols used here have already been defined. 

Per smooth spatial power density the performance index nay be written 
as 


T a Px2 - 3 


P.I. = / / [ ^(ff^+cyjg 4>- ff + a 2 (^ (a^OgF^ ~ |)) 2 ] dx dt 


o o 


She trial function for I is given as 


... 6.5 


*-*<,+ l l V cos 

it=1 ihp=1 


m 

2a 


6.6 


She trial function for Ng may be taken as 

L M K ^ trx m 

JSL = c V y b cos — * 

9 th ni, ® 2a 


6.7 


where c is a constant number, c can be determined from small tine 
step calculation and is assumed to be constant for the normal 
reactor operation time. 

Ijy using 6.2 and 6.3 in 6.5, the performance index changes to, 
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Using 6.9 for <J> in the performance index given by 6.8 it is found 
that 


I.X.-KV) *=’’ 2 1 

n =3 1 ^ 2 1 • • * jM 

1 .i 

low, by the help of constraint relation 6.4, the constrained coeffieiu. ts 

(b s n » 1,2,...,M) nay be substituted in terns of the independent 
v 1m 

coefficients (b ; £ = 2,3».».>Ii = 1,2,..,,M). Thus, perforr.ck.c- 
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index is a function of independent variables, 

* =2 ’ 3 1 

m = 1,2,. . .,M 


She coefficients can be determined by solving a system of 
nonlinear algebraic equations resulting from 




, %> ~ 0 & — 2, 3 >•••»! (C.io) 

m = 1,2,.».,M 

Onc e b ^ are solved, the rest part of the solution is quite easy. 

For the multiregions problems, the interface boundary conditions 
for the neutron flux and current are to be satisfied. These nay bo 
incorporated in the performance index of the problem in the following 
way : 


For a two region one dimensional problem let the fluxes are 
denoted by <^(x,t) and $ 2 (x,t). At the interface b 


^(b,t) = * 2 (b,t) 

3<J) (x,t) 

r\ J ■ 


1 3x 


= D 


3$ 2 (x,t) 
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6.11 


6.12 


x=b 


D ^ I> 2 being diffusion coefficients in the two regions. 

These equations 6.11 and 6.12 can be taken care of by includA v 
the terms 
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to the performance index. 

The reflector may be considered as an additional region and 
the above mentioned technique msy be used or an unreflected core with 
reflector savings may be treated. 

The methods adopted here may be extended for three dimensional 
physical reactor probl ems by assuming three dimensional trial functions 
for the state variables. Numerical evaluation of triple integrals 
will be involved. 
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